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We analyze the possibility that neutrino telescopes may provide an experimental determination 
of the slope A of the gluon distribution in the proton at momentum fractions x smaller than the 
accelerator reach. The method is based on a linear relation between A and the spectral index (slope) 
of the down-going atmospheric muon flux above 100 TeV, for which there is no background. Con- 
sidering the uncertainties in the charm production cross section and in the cosmic ray composition, 
we estimate the error on the measurement of A through this method, excluding the experimental 
error of the telescopes, to be ± 0.2. 



Atmospheric neutrinos and muons are the most important source of background for present and future neutrino 
telescopes, which are expected to open a new window in astronomy by detecting neutrinos from astrophysical sources. 

At energies above 1 TeV, atmospheric lepton fluxes have a prompt component consisting of neutrinos and muons 
created in semileptonic decays of charmed particles, as opposed to the conventional leptons coming from decays of 
pions and kaons. Thus a model for charm production and decays in the atmosphere is required. 

We base our model on QCD, the theoretically preferred model, to compute the charm production. We use a next- 
to-leading order perturbative QCD (NLO pQCD) calculation of charm production in the atmosphere, followed by a 
full simulation of particle cascades generated with PYTHIA routines |Q . 

We have already examined the prompt muon and neutrino fluxes in two previous papers (^|| (called GGV1 and 
GGV2 from now on). 

In our hrst paper ||, we found that the NLO pQCD approach produces fluxes in the bulk of older predictions (not 
based on pQCD) as well as of the recent pQCD semianalytical analysis of Pasquali, Reno and Sarcevic 0]. We also 
explained the reason of the low fluxes of the TIG model 0, the first to use pQCD in this context, which were due to 
the chosen extrapolation of the gluon partonic distribution function (PDF) at small momentum fractions x, and we 
confirmed the overall validity of their pioneering approach to the problem. 

In our second paper || , we analyzed in detail the dependence of the fluxes on the extrapolation of the gluon PDF 
at small x, which, according to theoretical models, is assumed to be a power law with exponent A, 



with A in the range 0-0.5. Particle physics experiments are yet unable to determine the value of A at x < 10 5 . We 
found that the choice of different values of A at x < 10~ 5 leads to a wide range of final background fluxes at energies 
above 10 5 GeV. 

Due to this result, in GGV2 we suggested the possibility of measuring A through the atmospheric leptonic fluxes 
at energies above 10 5 GeV, not the absolute fluxes, because of their large theoretical error, but rather through their 
spectral index (i.e. the "slope"). In particular, we now propose to use the slope of the flux of down-going muons. 

We want to stress that we are proposing to use down-going muons, at energies > 100 TeV, where prompt muons 
dominate over conventional ones, and not up-going neutrino- induced muons whose flux is orders of magnitude smaller. 
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While an important contribution to up-going muons is expected from astrophysical neutrinos, there is no background 
for down-going atmospheric muons. 

In this paper we further investigate the possibility of measuring A, in the more general context of an overall error 
analysis of our model. 

We can identify five potential causes of uncertainty in our final results. The first one is the possible presence of 
large logarithms of the type a s \np^ and a s lns (the latter are the so called "ln(l/x)" terms). The second is the 
treatment of the multiplicity in the production of cc at high energies. 

The third one consists of all the sources of uncertainty hidden in the treatment of particle cascades generated by 
PYTHIA. The fourth one is the uncertainty in the NLO pQCD charm production model we use. This includes the 
dependence of the fluxes on the three parameters of the model and the PDF's used. The fifth and final one is the 
choice of the primary cosmic ray flux, which is the input of our simulation. Of all these potential sources of uncertainty 
we conclude that only the last two are relevant. 

We deal with these five potential sources of error in turn. In Sect, [n] we address the question of the large logarithms 
a s lnpy and a s Ins, and in Sect. [II we analyze the problem of multiplicity in our charm production mechanism. 

In Sect. we consider the uncertainties due to the cascade generation by PYTHIA and to our NLO pQCD charm 
production (the core of our analysis). Here we evaluate the errors due to the parameters of the model, errors that 
affect the charm production cross section, the final differential and integrated fluxes and their spectral indices. We also 
determine how the final results (fluxes and their spectral indices) are affected by the choice of different extrapolations 
of the PDF's at x < lfT 5 . 

We are finally ready in Sect, [v] to discuss how A could be measured. We study the dependence on the different 
extrapolations of A at x < 10~ 5 , we consider the spectral indices and, using the discussion of Sect. IV, we provide an 
estimate of the errors on these indices and examine the feasibility of an experimental determination of A at x < 10~ a 
with neutrino telescopes. 

Finally in Sect. VI we discuss the error on the determination of A coming from the uncertainties in the elemental 
composition of the cosmic ray flux. 

The determination of A at small x < 10~ 5 is important because in this range saturation, unitarity and shadowing 
effects should become important. The PDF sets we use have been derived without including saturation effects. Even 
if this procedure seems to work very well in the HERA regime (where there might be some indications of saturation 
already, see e.g. ||), here we are extrapolating the resulting gluon PDF's at even smaller x values where saturation 
may become important. With respect to unitarity, using the expression of the Froissart bound on the gluon structure 
functions given in Eq. 31 of Ref. § we see that the extrapolated gluon PDF's we use, with A = 0.4 — 0.5, violate 
this bound at x values between 0.5 x 10~ 7 and 1 x 10~ 7 , for the characteristic momenta Q 2 ~ m 2 , ~ 3 GeV 2 we have, 
which corresponds to leptonic energies of 1 — 2 x 10 6 GeV. Always using the Froissart bound on the gluon PDF as 
given in the Eq. 31 of Ref. @, the gluon PDF's extrapolated with A < 0.3 encounter this bound at x < 10 -8 , what 
corresponds to leptonic energies larger than 10 s GeV, beyond the energy range relevant in this paper. Shadowing of 
the gluons in the atmospheric nucleons and nuclei, which we have not included here, could decrease the amplitude of 
the gluon PDF's by about 10% at x ~ 10~ 5 and up to as much as 30 to 40% at x ~ 10~ 7 (see e.g. ||), what would 
also change the effective value of A. There are no shadowing effects in the cosmic ray nucleons and nuclei. Only the 
dominant x of the gluons in the atmosphere is small in our calculation, while the dominant x of the partons in the 
cosmic rays is large. Thus shadowing effects do not depend on the unknown composition of cosmic rays, but could 
only be important for the atmospheric partons. The reason for the different characteristic values of x in the target 
and projectile partons is the following (for more details see GGV2). Due to the steep decrease with increasing energy 
of the incoming flux of cosmic rays, only the most energetic charm quarks produced count and those come from the 
interaction of projectile partons carrying a large fraction of the incoming nucleon momentum. Thus, the characteristic 
x of the projectile parton, x±, is large. It is x\ ~ O(10 -1 ). We can, then, inmediately understand that very small 
parton momentum fractions are needed in our calculation, because typical partonic center of mass energies y§ are 
close to the cc threshold, 2m c ~ 2 GeV, (since the differential cc production cross section decreases with increasing 
s) while the total center of mass energy squared is s = 2mi\rE (with mjy the nucleon mass, mjy — 1 GeV). Calling xi 
the momentum fraction of the target parton (in a nucleous of the atmosphere), then, x\x 2 = s/s — Am 2 /(2mNE) ~ 
GeV /E. Thus, X2 — O(GeV/0.1 E), where E is the energy per nucleon of the incoming cosmic ray in the lab. frame. 
The characteristic energy E c of the charm quark and the dominant leptonic energy Ei in the fluxes are Ei ~ E c ~ 0.1-E, 
thus x 2 ~ 0(GeV/ Ei). 
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II. IMPORTANCE OF THE a s ln(l/X) TERMS 



We address here a concern that has been expressed to us several times, about the applicability of perturbative QCD 
calculations, mostly done for accelerator physics, to the different kinematic domain of cosmic rays. 

Contrary to the case in accelerators, we do not have the uncertainty present in the differential cross sections || when 
the transverse momentum px is much larger than m c , uncertainty which is due to the presence of large logarithms of 
(p T + m^/m 2 .. The reason is that we do not have a forward cut in acceptance, and so the characteristic transverse 
charm momentum in our simulations is of the order of the charm mass, px — 0{m c ), and not px 3> 0(m c ) as in 
accelerator experiments. 

We may however, depending on the steepness of the gluon structure function A, have large logarithms of the type 
a s lns, known as "ln(l/x)" terms (here x ~ y/Arri^/s is the average value of the hadron energy fraction needed to 
produce the cc pair at hadronic center of mass energy squared s). These "ln(l/a;)" terms arise when the t-channel 
gluon exchange becomes large, and eventually they have to be resummed. Although techniques exist for resumming 



these logarithms 10 1, we have not done it. On the other hand we have phenomenologically altered the behavior of 
the parton distribution functions at small x by imposing a power law dependence of the form xf(x) ~ x~ x . This 
is analogous to resumming the ln(l/x) terms in a universal fashion and absorbing them in an improved evolution 
equation for the gluon density (such as the Balitskyii-Fadin-Kuraev-Lipatov (BFKL) evolution equation) , a 
procedure which increases A. For sufficiently large A, the large ln(l/x) terms should not be present. 

To find if our NLO cc cross sections are dominated by the ln(l/x) terms, we have used the following qualitative 
criterion flr2| . We have plotted the ratio 

R = aNLO ~ aLO (2) 
a LO a s \n(s/m 2 )/n 

as a function of the beam energy E. If the ratio is constant we are dominated by the ln(l/a;) terms and if it decreases 
we are not. The good behavior is a decreasing R. Figure 1 shows indeed that up to highest energy we consider in this 
paper, i.e. 10 11 GeV, R decreases for A > 0.2, but is roughly constant for smaller A's. This indicates that we are not 
dominated by the ln(l/ir) terms provided A > 0.2. 

Clearly, this test involving the R ratio does not say anything about ln(l/x) higher order corrections. One can only 
argue that if the ln(l/x) terms are not dominant at NLO (for R decreasing with energy) the corresponding [ln(\/x)] n 
terms may also be non-dominant in higher order corrections. In any event, the data on charm production that could 
be inferred at x < 10~ J , from the slope of atmospheric muon fluxes, really give information on the product of the 
gluon PDF and the parton cross section and a measurement of this product is useful. One can expect that the ln(l/x) 
terms at higher order may be better understood by the time the data will come. 



III. MULTIPLICITY IN CHARM PRODUCTION 



Another concern is the fact that at high energies the charm production cross section we use is sometimes larger 
than the total pp cross section. At first sight this seems absurd, but we show here that the cross section we use is the 
inclusive cross section, which contains the charm multiplicity, i.e. it counts the number of cc pairs produced, and so 
can be larger than the total cross section. On the other hand, the contribution of cc producing events to the total pp 
cross section, i.e. the cross section for producing at least one cc pair, is always smaller than the total pp cross section. 

We call ctqcd the perturbative QCD cross section of cc pair production in pp collisions, 

vqcd = ^2 ctqcd (ij -> cc) , (3) 

y 

where the sum is over the partons i and j in the colliding nucleons, and 

f dcr (ij — ^ cc\ 
VQCD(ij ->• cc) = J dx 1 dx 2 dQ 2 fi{ x l, Mjr)/j(lC2, Mf)- ( 4 ) 

Here da(ij — > cc)/dQ 2 is the ij — > cc parton scattering cross section, Q 2 is the four- momentum transfer squared, x t is 
the fraction of the momentum of the parent nucleon carried by parton i, and fi(x, fip) is the usual parton distribution 
function for parton momentum fraction x and factorization scale [ip. 

In the scattering of each pair of partons (one parton from the target and one from the projectile) only one cc pair 
may be produced, but the number of parton pairs interacting in each nucleon-nucleon collision is in general not limited 
to one and it increases with the number of partons f(x, n F )dx in each nucleon. 
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For A close to 0.5, cjqcd becomes larger than the total pp cross section a pv ~ 200 mb at E p ~ 10 10 GeV. It is 
obvious therefore that our results at high energy and large A are unphysical, unless multiplicity is taken into account. 
In fact, multiparton interactions should be taken into account already at a smaller cross section of order 10 mb, as 
determined in studies of double parton scattering Jl^ j . 

In order to incorporate multiparton scatterings into our analysis, we use an impact-parameter representation for the 
scattering amplitude, and ignore spin-dependent effects (cfr. |14[| ). Assuming the validity of factorization theorems, 
the mean number of parton-parton interactions ij — > cc in the collision of two protons at impact parameter b is given 
by 

n c5 (b) = £ / d 2 b'd Xl dx 2 dQ 2 da{l3 d Z ° C) Mf, b')f j (x 2 , b + V), (5) 

ij 

where fi(x, p, 2 ?, b)d 2 b is the number of partons i in the interval (x, x + dx) and in the transverse area element d 2 b at a 
distance b from the center of the proton. For simplicity of notation we drop the vector symbol in b and write b from 
now on. 

If n cS (b) -C 1, n cS (b) is the probability of producing a cc pair in a pp collision at impact parameter b. If n cS (b) > 1, 
n cS (b) is just the mean value of k, the number of cc pairs produced, at impact parameter b. Let the probability of k 
scatterings ij — > cc in a pp collision at impact parameter b be Pkcc{b). Then 

oo 

n cB (b) = J2 kp kcc{b)- (6) 

fc=0 

The /c-tuple parton cross section is obtained by integrating the probability of exactly k parton scatterings Pkcc(b) 
over the impact parameter b, 

the inclusive cross section for charm production is, thus, <7 c gmci = kokcc and the contribution of charm producing 
processes to the total cross section is a cS — ^2 k (Jkcc for k =/= 0. 

In our evaluation of charm production by cosmic ray interactions in the atmosphere, we must count the number of 
cc pairs produced in the pp collision. So we are precisely interested in the inclusive cross section <7 c gmci, which includes 
the number k of cc pairs produced per collision (the multiplicity). We hnd 



k " k 



:-t, , m ,, = kakcc = \ d 2 b22kP kc c(b) = I d 2 b n c5 {b). (8) 



This cross section can be larger than the total pp cross section, because it accounts for multiparton interactions. In 
particular, using cr cE , the contribution of charm producing processes to the total cross section defined above, the ratio 
Cccmci/ccc gives the average charm multiplicity. 

Notice that here we consider only independent production of cc pairs, so that from each pair of colliding partons it 
results only one cc pair, and we neglect coherent production of multiple cc pairs in 2— >4, 2— >6, etc. processes. This 
will underestimate the charm production cross section. 

We assume in the following that the partonic distributions fi(x, Hp,b) factorize as 

fi(x,fi 2 F ,b) = fi(x,fi%)pi(b), (9) 

where fi(x,[i 2 F ) is the usual parton distribution function, and p%(b) is the probability density of finding a parton in 
the area d 2 b at impact parameter b. We normalize Pi(b) to J d 2 b Pi(b) — 1, to maintain the usual normalization 
/ dx xfi (x) = 1 . The factorization in Eq. (j^) is consistent with the usual parton picture and with our assumption of 
no parton-parton correlations. 

The mean number of ij — > cc scatterings at impact parameter b then becomes 

n c s(b) = ^2atj(b)<J Q cD(ij -> cc), (10) 

ij 

where 
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aij (b) = / dH'pi^pjib + b 1 ) 



(ii) 



is an overlap integral, and (JQCD^j - > cc) is the QCD parton-parton cross section for ij — » cc, as in Eq. (Q). From 
the normalization of pi(6) it follows that J (i 2 6 (%■(&) = 1 for every Hence from Eqs. (§) and © we find 

Cccincl = &QCD, (12) 

where uqcd, given is Eq. (||), is the charm production cross section calculated within QCD. This justifies our use of 
cqcd a $ Cccinci in the calculation of the atmospheric fluxes. 

The way in which we use cr C cmci in our simulation is as follows. We input only one cc pair per pp collision at a given 
energy E, and multiply the outcome by cr c emci, which includes the cc multiplicity. We make, therefore, the following 
approximation in the kinematics of the cc pairs produced in the same pp interaction. Even if in a real multiparton 
collision the energy available to the second and other cc pairs is smaller than E, we are neglecting this difference. 
This is a very good approximation because the fraction of center of mass energy that goes into a cc pair is of the 
order of y/IJs ~ y/10 GeV/E < 1 at the high energies we are concerned with. 

Although we have explicitly proven Eq. (|l2|) in the absence of parton-parton correlations, the same result can be 
obtained when correlations are present (see sect. V of Ref. fj"5f| and references therein). What is proven even in the 
presence of correlations is that the pQCD single scattering cross section oqcd is equal to the average number of 
parton-parton collisions, call it < N >, multiplied by the contribution of cc producing events to the total cross section 
(the cross-section cr c5 defined above), namely ctqcd =< N > a cS (while the inclusive cross section is equal to the 
average multiplicity of cc pairs multiplied by <r ce ). < N > may in general contain contributions from two types of 
collisions. One type consists of collisions of pairs of partons (consisting of one parton from each interacting nucleon) 
which interact only once at different points of the transverse plane. Each collision of this type results in our case in 
one cc pair produced. The second type consists of rescatterings in which one parton of one of the nucleons and its 
interaction products interact with several partons of the other nucleon. In interactions of this type, which are much 
rarer than the first ones, the number of cc pairs produced not necessarily coincides with the number of collisions. If 
rescatterings can be neglected, then < N > is the average number of cc pairs produced and ctqcd is the inclusive 
cc production cross section as stated in Eq. (|l2|). Otherwise small rescattering corrections, to our knowledge not yet 
calculated |]l6| , are necessary (rescatterings would also modify the energy spectrum of the particles produced). 



IV. UNCERTAINTIES DUE TO CASCADE SIMULATION, PARAMETERS OF CHARM PRODUCTION 

MODEL AND CHOICE OF PDF'S 

In our first paper (GGV1) we considered the uncertainties related to the cascade generation in PYTHIA. There 
we tried different modes of cascade generation, different options allowed by PYTHIA in the various stages of parton 
showering, hadronization, interactions and decays, etc., without noticing substantial changes in the final results 
(differing at most by 10 %). These uncertainties are however very difficult to quantify, due to the nature of the 
PYTHIA routines. Since these uncertainties are small, we neglect them in this analysis and continue to use PYTHIA 
with the options described in GGV1 as our main choice for the simulation: 'single' mode with showering, 'independent' 
fragmentation, interactions and semileptonic decays according to TIG. 

In our 'single' mode we enter only one c quark in the particle list of PYTHIA, and we multiply the result by a 
factor of 2 to account for the initial c quark. PYTHIA performs the showering, standard independent fragmentation, 
and follows all the interactions and decays using default parameters and options. In GGV1 we have argued that this 
'single' approach is equivalent to what we called 'double' mode, in which both cc partons are placed in the initial 
event list, in the first step of a standard cascade evolution. The 'single' option is chosen thus, because it reduces 
considerably the computing time. 

Important sources of uncertainty are contained in our charm production model, which is NLO pQCD as implemented 
in the MNR program ||, calibrated at low energies. 

The calibration procedure consisted in the following: 

• choosing a PDF set from those available and fixing the related value of Aqcd ;Q 



1 We note that Aqcd can be chosen in the MNR program independently of the PDF and therefore can constitute an additional 
independent parameter of our model. We have opted however to choose the value of Aqcd assumed in the PDF set being used, 
for consistency. 
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• choosing m c , fip and iir, which are the charm mass, the factorization scale and the renormalization scale 
respectively, so as to fit simultaneously both the total and differential cross sections from existing fixed target charm 
production experiments Jl7],[l8| at the energy of 250 GeV, without additional normalization factors; 

• checking that the total cross section generated after the previous choices fits reasonably well the other existing 



experimental points for fixed target charm production experiments 19 1. 

Besides the choice of the PDF set, our procedure has the freedom to choose reasonable values of the three parameters 
m c , [If, and iir so as to fit the experimental data. In GGV1 and GGV2 we made the standard choice pjl9|] of 



[i F = 2m T , n R = m T , (13) 

where mj = \Jp\ + rric is the transverse mass. The values of the charm mass are taken slightly different for each 
PDF set, namely: 

m c = 1.185 GeV for MRS Rl, (14) 

m c = 1.310 GeV for MRS R2, (15) 

m c = 1.270 GeV for CTEQ 4M, (16) 

m c = 1.250 GeV for MRST. (17) 

Here we explore the changes induced in cross sections and fluxes at high energies by different choices of m c , ^f, 
and /in which fulfil our calibration requirements. 



We have performed this analysis with the most recent PDF set: the MRST |20| (other PDF's give similar results). 

ibration procedure described above 



At first we fix A = and then we examine other values of A. We note that the ca 
is independent of A because it involves only relatively low energies, where the low x extrapolation is not an issue 



A. MRST A = 0: fluxes 



We considered the A = case because it is the most significant for the evaluation of the uncertainties in the spectral 
indices, as it will be clear in the next subsection. We have considered the following reasonable ranges of the parameters 

1.1 GeV < m c < 1.4 GeV, (18) 
0.5 m T < ii F < 2 m T , (19) 
0.5 to t < fi R < 2 tot, (20) 



where the bounds on m c come from the 1998 Review of Particle Physics 21 1, while those for tip and i±r are those 
used in the existing literature J|[l5| . 

Within these ranges we have looked for values of the three parameters capable of reproducing the experimental data 
in our calibration procedure described before. Table [| summarizes the different sets of parameters: we have varied 
the charm mass through the values m c — 1.1, 1.2, 1.25, 1.3, 1.4 GeV (m c = 1-25 GeV was our previous optimal choice 
for MRST in Eq. (|l7|)) and then, for each value of m c , we have found values of the factorization and renormalization 
scales that reproduce the experimental total cross section cr ce = 13.5 ± 2.2 iih at 250 GeV JIt]]. In particular, for each 
value of m c , we took Lip = tot/2, tot, 2tot and then, for each m c , tMp pair, found the value of iir which best fits the 
data (see Table 0). 

We have also checked that these choices give good fits to the differential, besides the total, cross sections at 250 GeV 
without additional normalization factors, as done for the original choice of parameters in GGV1. For m c = 1.1 
GeV we had to choose values of ixr slightly outside the range in Eq. ( |20| ) (but we have kept these values in our analysis 
anyway) . 

For all the sets of parameters in Table Qwe have run our full simulations for the MRST, A = case and the results 
are described in Figs. 2-4. 

In Figs. 2a and 2b we show the resulting total charm production cross section a c5 for all of the fifteen sets of 
parameters in Table |[ together with recent experimental data (from Table 1 of Ref. JO]] , where all the data for pp or 
pN collisions have been transformed into a a cS cross section following the procedure described in GGV1). Fig. 2b is 
an enlargement of the region of Fig. 2a containing the experimental data. 

In Fig. 2a we see the spread of the cross sections, which is more than one order of magnitude at 10 11 GeV. Above 
250 GeV, one can clearly distinguish three "bands" of increasing cross sections for fip = tot/2, tot and 2tot- Within 
each "band" the cross sections increase with increasing values of m c (and correspondingly smaller values of (J-r), in 
Table Q. Our standard choice (m c — 1.25 GeV, fip = 2tot, [ir = rap) proves to be one of the highest cross sections 
we obtain. 
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In Fig. 2b we see better how all of these cross sections verify our calibration procedure. They pass through the 
point at 250 GeV fL7j , agree with the point at 400 GeV J2^] and disagree with the very low experimental point at 200 
GeV |2^]. The lower values of fip — my/2 and mj fit better the lowest experimental point at 800 GeV p3], while 
the higher value of \xf — 2rriT fits better the upper point at 800 GeV j2q| . 

We believe that the spread of the total cross sections shown in Fig. 2a provides a reasonable estimate of the 
uncertainty of our charm production model at fixed A. Since our standard choice of parameters (m c = 1.25 GeV, 
fiF — 2my, [ir — tut) gives one of the highest cross sections (in better agreement with the more recent value of the 
cross section at 800 GeV ) , the uncertainty band should be added under each of the cross section curves calculated 
with our standard choice of parameters (like the curves shown in Fig. 1 of GGV2). 

Fig. 3 illustrates the corresponding spread of the final prompt fluxes. Although our results are for the MRST PDF's 
extrapolated with A = (the value of A which gives the lowest fluxes) similar spreads result from other PDF's and 
A's. We show the flux of muons; the fluxes of muon-neutrinos and electron-neutrinos are essentially the same. 

Similarly to what happens with cross sections in Fig. 2, the fluxes in Fig. 3 increase with fJ-F = wt/2, tot and 
2tot and, within each band, they increase with increasing m c (and correspondingly smaller values of hr), in Table |. 
At energies around 10 6 GeV the total uncertainty is almost one order of magnitude and decreases slightly for higher 
energies. If we would decide to work only with [ip = 2tot (which fits the experimental measurement at 800 GeV 
with the highest cross section) , the uncertainty would be greatly reduced: the fluxes in this rather narrow band differ 
by less than 40%. We observe that the flux calculated with our standard choice of parameters (m c = 1-25 GeV, 
fip — 2to^, /ifl = rap) is almost the highest, as was the case for the corresponding cross section in Fig. 2. 

In Fig. 3 we also indicate the conventional and prompt fluxes from TIG; we notice that the TIG prompt flux is 
within our band of uncertainty, which is reasonable since TIG used a low A = 0.08 value for their predictions (see the 
discussions in GGV1 and GGV2). 

B. MRST A = 0: spectral index 

In our previous paper GGV2, we pointed out that an experimental measurement of the slope of the atmospheric 
lepton fluxes at energies where the prompt component dominates over the conventional one, might give information on 
the value of A, the slope of the gluon PDF at small x. The best flux for this measurement is that of down-going muons, 
because the prompt neutrinos have first to convert into muons or electrons through a charged current interaction in 
order to be detectable in a neutrino telescope. 

In this section and in the following two we consider the uncertainties in our method to determine A. In this section 
we examine those coming from the charm production model, in Sect. [v| those related to the non linearity of our model, 



and in Sect. VI those coming from the unknown composition of the cosmic rays at high energies. 

The slope of the fluxes or spectral index is ae(Ei) — —d\n(f)e(Ei)/dlnEi, with i = /i ± , + or v e + v e . In other 
words, the final lepton fluxes are 

UEi) oc £- a ' (Sl) (21) 
In GGV2 we found a simple linear dependence of ae on A, namely 

a e (E e ) = b e (E e , 7 , A) - A ~ b e {E e ) - A, (22) 

where b^iEf) is an energy dependent coefficient evaluated using our simulation with A = and fixed 7. As argued in 
GGV2 (cfr. Eqs. (35) and (36) therein), the coefficient bi(Ei,j, A) depends mildly on A and can be well approximated 
by its value for A = (see Sect. |y|). The coefficient be(Ei,j, A) depends almost linearly on 7, the spectral index of 
the primary cosmic rays. We recall that the equivalent nucleon flux for primary cosmic rays is expressed as 

<j) N {E) ozE-1- 1 . (23) 

The linear dependence of 6^(^,7, A) on 7 can be written as 

bi{Eg, 7, A) = b~i(Ei, 7, A) -I- 7, (24) 

where bg(Ei, 7, A) depends mildly on A and 7 as we will prove in Sect. [v| and Sect. VI, respectively. 



J We have included in be the +1 term coming from the —1 in the exponent of Eq. 
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Given bt(Ee) from our model, an experimental measurement of ae at energy Eg would immediately give A corre- 
sponding to a value of x ~ GeV /Eg, as we discussed in GGV2. A measurement at Eg ~ 10 6 GeV = 1 PeV would give 
A at x ~ 10 -6 , a value of x unattainable by present experiments. 

For the time being we keep fixed the value of 7 (7 = 1.7 below the knee, and 7 = 2.0 above the knee, as in GGV1 
and GGV2); only in Sect. [?] we will consider changes in the value of 7. 

The feasibility of a measurement of A depends, therefore, on the uncertainties in bp{Ep). Here we discuss those 
coming from the charm production model. 

Fig. 4 shows the —b p corresponding to the fluxes of Fig. 3 as functions of E^. In the region of interest E^ > 10 5 -10 6 
GeV, the values of — within each "band" decrease with increasing m c (and correspondingly smaller values of /j-r), 
in Table § 

The spread of b^ due to the choice of fj,p > Mi? an d fn c is Ab^ ~ 0.1, or Ab^/b^ ~ 0.03, much smaller than 
the uncertainty A^/^ ~ 10 of the absolute flux in Fig. 3. If we would for some reason restrict ourselves to the 
= 2 rriT band, the uncertainty on b^ would become even smaller, A6 M ~ 0.03. We will refer to this error as Ab par 
in the following, as it is related to the choice of parameters in the charm production model, and consider half of the 
spread in Fig. 4 to evaluate it. Therefore 

Ab par ~ 0.05 (0.015), (25) 
where the value in parenthesis corresponds to the /j,p = Ivtit band. 



C. MRST A = A(T) 

So far we used A = only. This case determines the uncertainty of the bt{Et) function which enters in the 
determination of A through the atmospheric muon fluxes. 

Here we study an "intermediate" value of A. We continue to use the MRST PDF, but with the value of A = A(T) 
given by the slope of the lowest tabulated value of x (see GGV2 for more explanations). This value depends on Q 2 
and is about 0.2-0.3. 



We repeat the same analysis of subsection IV A. However, for simplicity, we report the results for four selected sets 
of values for the parameters in Table |. The first set (m c — 1.1 GeV, hf = 0.5 my, [ir = 2.53 my) gives a lower 
bound for the fluxes. The second set [m c = 1.4 GeV, hf = 2m/r, M-R = 0-61 mr) gives an upper bound for the fluxes. 
The remaining two sets are cases in the fiF = ^mp "band". 

The results are plotted in Fig. 5. The general features of Fig. 5 coincide with those of Fig. 3, except for an overall 
increase in all the fluxes due to the larger value of A. The total spread of the fluxes given by the two limiting cases, 
as well as the spread within the narrower /j,p = 2rriT band, are comparable to those found for A = 0. As in Fig. 3, 
our standard choice of parameters (m c = 1.25 GeV, /Up = 2mp, (J>r = 1.0 my) yields almost the highest flux. 

We conclude that similar features would be obtained for other values of A: our "standard choice" flux would be 
almost the highest in a band of uncertainty whose width is similar for all values of A. The fluxes in the uncertainty 
band of Fig. 5 are consistent with older predictions (see GGV1 and references therein) and with the prediction by 
L. Pasquali et al. |Q. 



D. Other PDF's 



Another source of uncertainty for the final fluxes and spectral indices is the choice of the PDF set. As in GGV2, we 
consider here four recent sets of PDF's: MRS R1-R2 [g, CTEQ 4M (27) and MRST @, with the standard choice 
of parameters of Eqs. @,@,(||), (H),(0). 

Figs. 6a and 6b show the muon fluxes (top panels) and spectral indices (bottom panels) for the two limiting cases 
of A = (Fig. 6a) and A = 0.5 (Fig. 6b). In both cases the n fluxes show at most a 30 — 50% variation depending 
on the PDF used. The uncertainty in the spectral indices for E^ > 10 5 - 10 6 GeV is A6 M < 0.02, or Ab^/b^ < 0.01. 
This error will be denoted as AbpuF in the following, namely (again dividing the spread by 2) 

A&P£> F ~0.01. (26) 

These uncertainties, related to the PDF's, are smaller that those due to the choices of mass scales (see Figs. 3-4). 
We conclude that, provided different PDF's are calibrated in a similar way (i.e. same values of jip, fj,R and m c , chosen 
to fit the experimental data of our calibration), the final fluxes and spectral indices are very similar. The main source 
of uncertainty resides therefore in the choice of the mass parameters, rather than the adoption of a certain PDF set. 
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V. DETERMINATION OF A WITH NEUTRINO TELESCOPES 



In GGV2 we have given a detailed analysis of the dependence of the final fluxes and spectral indices on A for 
different PDF's. In this section we show that the spread in our results due to A is larger than the one due to the 
choice of fj,F, Hn., mc and of the PDF set, analyzed in the previous section. This is good news for the possibility of 
measuring A, since the spread in a M , due to different A's, is the signal we want to detect, while the spread due to 
other factors constitutes the theoretical error of this measurement. 

Figs. 7-10 show how the \i flux and its spectral index depend on A. We used MRST with variable A = 
0, 0.1, 0.2, 0.3, 0.4, 0.5 and our standard choice of parameters (m c = 1.25 GeV, hf = 2mr, fiR — 1.0 mj). 

Fig. 7 contains the differential muon fluxes. At the highest energies the fi fluxes are spread over almost two orders 
of magnitude. To each of the curves in this plot we need to assign a band of uncertainty of about one order of 
magnitude coming from the choice of the PDF and of the parameters m c , /if, and fiR (see Fig. 3). Thus the curves 
become entirely superposed with each other. This makes it impossible to derive the value of A from an experimental 
measurement of the absolute level of the fluxes. However, the uncertainties in the spectral index of these prompt 
muons are much smaller and a determination of A becomes possible using the slope of the muon fluxes instead of their 
absolute level. 

Fig. 8 shows the -E 2 -weighted integrated fluxes as functions of the muon energy. The slant lines indicate the number 
of particles traversing a km 3 detector over a 2ir sr solid angle. Even for the highest predicted fluxes, less than 1 particle 
per year will traverse the km 3 detector for energies above 10 s GeV. Moreover, while prompt muons can be detected 
directly, prompt neutrinos have first to convert into charged leptons before being detected. The smallness of the 
charged current interaction effecting the conversion considerably lowers the detection rate of neutrinos. Therefore, 
the slope of the charm component of the atmospheric leptons can be studied in neutrino telescopes only using 
atmospheric muons coming from above the horizon, and only in a narrow range of energies, between a lower limit of 
~ 10 5 — 10 6 GeV, above which the prompt component dominates over the conventional one, and an upper limit 
of ~ 10 7 — 10 s GeV, above which the detection rates become negligible. 

In practice, the spectral index of the prompt muon flux may be estimated by the difference of two integrated muon 
fluxes above two different energies, e.g. 10 6 and 10 7 GeV. 

Figs. 9, 10 prove the validity in our model of Eq. (p2f), which is at{Et) = bi(Ei) — A. In Fig. 9 we plot the 
spectral indices —ai(Eg) for the different values of A, both as directly calculated with our simulation (solid lines) and 
as —bi(Ei) + A (dotted lines), where be(Ee) is ag with A = 0. The two almost coincide, in the interval of interest, 
Et > 10 6 GeV. Their difference, a e {E e ) - b e (E e ) + A, given in Fig. 10, is small, about ~ 0.03 at E > 10 6 GeV. 
This difference stems from the mild dependence of bi(Eg) on A and need to be added to the the other uncertainties 



evaluated in Sect. |Vj. We will refer to this error, due to the non linearity in A of Eq. (|22|), as 

~ 0.015, (27) 

which again is evaluated dividing by 2 the spread in Fig. 10. 

We see in Fig. 9 that AA ~ Aa, therefore we would need an uncertainty in the spectral index of order 0.1 to 
determine A with the same accuracy. We will show now that this is roughly the uncertainty related to our theoretical 
model. 

In fact, we can finally estimate the total uncertainty in the determination of A coming from our theoretical model 
(that is, excluding the uncertainty due to the unknown composition of cosmic rays). We sum together the three 
spreads of bi(Ei) previously calculated in Eqs. (p5h, ( p6| ) and (p7[), to obtain the final uncertainty from the charm 
production model, 

{AX) charm ~ (A6 M ) cterm ~ 0.075 (0.04), (28) 

where the number in parenthesis corresponds to fixing fip — 2m,T in the charm model. 

If the theoretical uncertainties so far presented would be the only ones affecting the determination of A through a 
measurement of the slope of the down-going muon flux, we could expect to get to know A with an uncertainty of about 
A A ~ 0.1. However, even excluding experimental uncertainties in the neutrino telescopes themselves, the uncertainty 
increases when our ignorance of the composition of the cosmic rays at high energy is taken into account, as we show 
in the following section. 



5 We summed the errors linearly. Summing in quadrature would give (AX) c h a rm — (A& M ) c f l(lrm ~ 0.053 (0.023). 
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VI. UNCERTAINTY FROM COSMIC RAY COMPOSITION 



The final uncertainty we consider in the determination of A comes from the poorly known elemental composition 
of the high energy cosmic rays. 

The spectral index of the cosmic rays 7 enters almost linearly in the slope of the atmospheric leptons. From 
Eqs. (H) and (§|) we have 

oti{Et) = b t (Eg, 7, A) + 7 - A. (29) 



So far we have kept 7 fixed, thus the uncertainty A& M calculated in Eq. (|28|) is actually an uncertainty in b^. We are 
going now to evaluate the uncertainty due to 7. 

The non-linearity of Eq. (|2^) with respect to 7 is mild, as we have argued analytically in Sect. V of GGV2 and 
we show here using our numerical simulation. We have conducted a fe w tri al runs of our simulation simply changing 



the values of 7 used for the primary flux. We recall from subsection |IVE| that in our model we used 7 = 1.7, 2.0, 
respectively below and above the knee at E = 5 10 6 GeV . We have run our simulation changing these values of 7 by 
±0.1, ± 0.2 ^, both above and below the knee, to see the error produced when taking bg computed at fixed 7 (our 
usual values) in Eq. (^) and thus leaving a pure linear dependence on 7. We used the MRST PDF, with A — 0, but 
similar results are obtained with other PDF's and A's. 

In Fig. 11a we plot the spectral index —ctg(Eg) for the different values of 7, both as directly calculated with our 
simulation (solid lines) and as —bg(Eg; 7 = 1.7, 2.0; A = 0) — 7 (dotted lines), i.e. using our standard values for 7 of 
the primary flux and adding an increment in 7 equal to ±0.1, ± 0.2. In this way the "central value" of these spectral 
indices corresponds to the A = case of Fig. 9. We can see that the dotted and the solid lines almost coincide, 
especially in the interval of interest for Eg > 10 5 — 10 6 GeV, proving the validity of Eq. (p9[). The uncertainty in be 
due to this non-linearity, that we call A7 nors _;i n , evaluated in terms of the difference ctg — bg — 7, is plotted in Fig. lib 
and, in the energy range of interest, is 

A7„ on _ Mn ~ 0.02. (30) 

We will now consider the error due to the poorly known elemental composition of the high energy cosmic rays. 
Concerning charm production, the relevant cosmic ray flux to be considered is the equivalent flux of nuclcons impinging 
on the atmosphere. For a given cosmic ray flux, the equivalent flux of nucleons 4> eq (E^) depends in general on the 
composition of the cosmic rays. Nuclei with atomic number A and energy Ea, coming with a flux (Pa{Ea), contribute 
an amount A4>a{AEn) to the equivalent flux of nucleons at energy En = Ea/A. So in total 



->eq 



{E N ) = J2 A MAE N ). (31) 



The uncertainties in the equivalent nucleon flux arise from the poorly known composition of cosmic rays in the energy 
range above the so-called knee, Ea ~ 10 6 GeV. 

The actual 7 that enters into our proposed determination of A is the spectral index of the equivalent nucleon flux 
7e g , the equivalent cosmic rays spectral index for short. The equivalent nucleon flux is written as 4> eq oc E N 9 , so 
that the spectral index ^ eq is given by 

7e 9 + l = -f^ = -^£ AMlA + 1), (32) 
4> eq 3E N 4> eq 

where ^a is the spectral index of the component of atomic number A, i.e. 4>a{Ea) — k,AE~^ <A . 

We have calculated <j> eq and j eq using the experimental data of JACEE CAS A- MIA @, HEGRA [|o(, and 
the data collected by Biermann et al, in Table 1 of Ref. [[Jl], each with their respective compositions. Figs. 12 and 13 
show the 4> eq and the j eq so obtained. Only the data of CASA-MIA |2jJ and HEGRA Q reach energies E N < 10 8 
GeV, so we have not extended our analysis beyond 10 8 GeV. 

We have calculated the error band associated to 7^ in two different ways, because of the different parametrization 



of the composition used in Refs. [£8| to |31|. Refs. p8pl[ give separate power law fits to the spectrum of each cosmic 
ray component, 



'Notice that these values of 7 are some of the most probable values (see Fig. 13). 
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J a (33) 
where the parameters Ua and 7a have errors AkA and A 7a- Standard propagation of errors gives, in this case, 



and 



A 7e 



E 



A(f> e 



<f>% 



E^ 



AkA 

k A 



(ln(AE N )A 7/ 



1/2 



(34) 



(lA ~leqf 



AkA 

k A 



1 - 



(lA ~ leg) ln(AE N 



1/2 



(35) 



where <f> A is evaluated at Ea — AEn- 

Refs. [p9pC|], give a power law fit to the total particle flux 



and a composition ratio ta(Ea) in terms of which 

MEa) -- 



kE A ^ 



r A (E A ) <f>(E A ). 



(36) 



(37) 



These experiments distinguish only between a light and a heavy component. We assign atomic number 1 to the 
light component and 56 to the heavy one (which we call "iron" ) . Here fc, 7, and r A have errors Ak, A7, and Ar Al 
respectively. The equivalent nucleon flux is still given by Eq. ( |3l| ) , while standard propagation of errors gives in this 
case 




\ 2 

Xj + \}l A( t > AMAE N )A lA 



1/2 



(38) 



We omit the much longer expression for A7 eg . For simplicity, we have neglected the error coming from the energy 
dependence of r A , which we expect to be much smaller than the others. In Fig. 12 we show the equivalent nucleon flux 
<j) eq . It is clear from the figure that the systematic uncertainties dominate, with spreads between different experiments 
of up to a factor of 4. 

The uncertainties in the equivalent spectral index "f eq are smaller, as can be seen in Fig. 13, where only HEGRA 
and CASA-MIA extend to the energy region above the knee which is important to us. 

We can consider, for example, an energy En — 10 7 GeV, which is likely to determine the leptonic fluxes at around 
Ei ~ 10 6 GeV , energy at which we would like to measure A through the spectral index (we recall from GGV2 that 
E e < 0.1 E N ). 

At this energy En, from Fig. 13, we may take half the difference between the central values of the CASA-MIA and 
HEGRA data as an indication of the systematic uncertainty on 7 eg , 



A 7 . 



syst 



0.1. 



(39) 



Using the CASA-MIA data and the related error band, instead of the very spread HEGRA data, we can expect a 
reasonable statistical uncertainty 



A 7stat ~ 0.05. 



(40) 



Since oti depends linearly on ^ eq and A, the same uncertainties apply to a determination of A. The total uncertainty 
in the determination of A coming from the unknown composition of cosmic rays is now simply the sum of Eqs. (j3C" 
© and (HT 



(AA) e 



(A 7 eg)c 



if summing the errors linearly, or 



(AA) 



comp 



(A 7 e 9 ) 



q J comp 



0.17, 



0.11, 



(41) 
(42) 
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if we sum them in quadrature. 

Finally, we can now combine all the uncertainties together, to compute the overall theoretical error in the determi- 
nation of A with neutrino telescopes. From Eqs. (§|), (§|), ©, @, @, and @ we obtain 

AA-0.25 (0.21) (43) 

if summing errors linearly, or 

AA-0.13 (0.11), (44) 
if summing in quadrature, where the numbers in parenthesis correspond to the fip — 2to^ "band" in the charm model. 



VII. CONCLUSIONS 



We have examined in detail the possibility of determining the slope A of the gluon PDF, at momentum fraction 
x < 10~ 5 , not reachable in laboratories, through the measurement in neutrino telescopes of the slope of down-going 
muon fluxes at ~ x^ 1 GeV . 

The method we are proposing may reasonably well reach x ~ 10~ 7 , what would require 10 PeV in muon energy. 
At this energy, there would still be 50 events from charm if A = 0.5 and 10 events if A = 0. But the best measurement 
could be done between 100 TeV and 1 PeV of muon energy, i.e. between x ~ 10~ 4 and x ~ 10~ 6 . Present data do 
not go below x ~ 10~ 5 and the Large Hadron Collider (LHC) will not do any better. The reason is that the dominant 
values of x in the production of a heavy particle of mass M and rapidity y are of order x ~ [Mexp (±y)/y/s] (see 
for example p2|) where yfs is the center-of-mass energy of the hadron collision. Thus the smaller values of x are 
obtained with the smaller M and larger y for fixed s/s (14 TeV at the LHC). Even if exhaustive studies of the possible 
minimum x to be reached at the LHC have not yet been carried out p3[ | , it is known that the experiments will explore 
the central rapidity region (the CMS and ATLAS detectors will cover y < 0.9 only) and that bottom can be tagged, 
but most likely not charm ^ |34| . This means that the lowest x that LHC is expected to reach, assuming realistically 
that charm will not be tagged, is x ~ mj exp (— 0.9)/t/s = 1.5 x 10 -4 . Therefore, the method proposed here may give 
information on the gluon PDF at x < 10 -5 , a range not reachable in laboratory experiments in the near future. 

To this end we studied the dependence of the leptonic fluxes and their slopes on A. The slopes depend almost 
linearly on A. We studied the uncertainties of the method we propose (excluding the experimental errors of the 
telescopes themselves). These come mainly from two sources: the free parameters of the NLO QCD calculation of 
charm production and the poorly known composition of cosmic rays at high energies. 

We have seen that, for a fixed value of A, the uncertainties give rise to an error band for the leptonic fluxes of 
almost one order of magnitude at the highest energies. This makes impossible a determination of A based solely on 
the absolute values of the fluxes, therefore we propose using the slopes of the fluxes. In particular we are proposing 
to use down-going muons, for energies > 100 TeV , where prompt muons dominate over conventional ones, and 
not up-going neutrino-induced muons whose flux is orders of magnitude smaller. While an important contribution to 
up-going muons is expected from astrophysical neutrinos, there is no background for down-going atmospheric muons. 

The overall theoretical error, from the charm production model, on the measurement of A, is (AA) c / larm < 0.10. A 
comparable error, due to uncertainties in the cosmic ray composition, (AA) comp < 0.15, must be added, so that the 
overall error in the determination of A with neutrino telescopes is AA ~ 0.2. 

These errors may be reduced by improving the experimental knowledge of the charm production cross sections and 
of the cosmic ray composition around and above the knee. 



5 Tagging is done by finding the point where the quark decays (called the vertex). The probability of decaying is exponential, 
and higher in the region close to the collision point. The only way to tell charm from bottom is by the distance from the 
collision to the vertex and bottom quark lives longer than charm. Thus, to detect charm with a good degree of confidence, one 
needs to select vertices close to the collision point. But in this region the vertices from bottom decay dominate, because the 
number of decay channels of the b quark is five times larger than that of the c quark. 
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FIGURE CAPTIONS 



Fig. 1 The ratio R = {<7nlo — &lo)/ (clo^s rn ( s / m c)/ 7r ) is plotted as a function of the beam energy E, for the 
different values of A used with the MRST PDF. 

Fig. 2 Total cross sections for charm production a C c, up to NLO, calculated with MRST (A = 0) and the values of 
m c , fip, fin of Table |[ are compared with recent experimental values ]l7| , [T9| , p2| -p5| . For each "band" in the 
figures (i.e. for fip = m-r/2, mr and Imp) the cross sections increase with increasing m c (and correspondingly 
smaller values of fin) in Table | (Fig. 2b is an enlargement of Fig. 2a). 

Fig. 3 Results for MRST A = 0. The i? 3 -weighted vertical prompt fluxes, at NLO, are calculated using the values 
of m c , fip, fin of Table ||and compared to the TIG B conventional and prompt fluxes. For each "band" in the 
figure (i.e. for ftp = tot/2, mp and Imp) the fluxes increase with increasing m c (and correspondingly smaller 
values of fi R ) in Table g. 

Fig. 4 Spectral indices — fe M of the fluxes plotted in Fig. 3, for the MRST A = case. For each "band" in the figure 
(i.e. for fip = mp/2, mj and Imp) the spectral indices decrease with increasing m c (and correspondingly 
smaller values of fin) in Table |. 

Fig. 5 Results for MRST A = A(T). The ^-weighted vertical prompt fluxes, at NLO, are calculated using selected 
values of m c , fip, fin from Table || and compared to the TIG Q conventional and prompt fluxes. 

Fig. 6 Results for MRS R1-R2, CTEQ 4M, MRST, for A = (Fig. 6a) and A = 0.5 (Fig. 6b), with standard choice of 
parameters m c , fip, fin- Top part: ^-weighted vertical prompt fluxes, at NLO. Bottom part: related spectral 
indices — a M (for the A = case, — a M = —b^). 

Fig. 7 Results for MRST A = — 0.5 (solid lines). The ^-weighted vertical prompt fluxes, at NLO, are compared 
to the TIG || conventional and prompt fluxes (dashed lines). 

Fig. 8 Results for MRST A = — 0.5 (solid lines). The ^-weighted integrated vertical prompt fluxes, at NLO, are 
compared to the number of particles traversing a km 3 2ir sr detector per year (dotted lines). 

Fig. 9 Results for MRST A = — 0.5. The spectral indices —ai(Ei) for the different values of A, calculated directly 
by our simulation (solid lines) are compared to the corresponding terms —bi{Ei) + A (dotted lines). 

Fig. 10 Results for MRST A = — 0.5 (solid lines). The error of Eq. (|2|) is evaluated in terms of the difference 
at(Et) - bt(Et) + \. 

Fig. 11 Results for MRST A = for different values of 7. a) The spectral indices —ag(Ei) for the different values 
of 7, calculated directly by our simulation (solid lines) are compared to the corresponding terms —bi(Et;^i = 
1.7, 2.0; A = 0) — 7, with increments in 7 equal to ±0.1, ± 0.2 (dotted lines). The curves arc labelled by the 
related value of 7 above the knee (7 = 2.0 is our "standard value"), b) Uncertainty due to the non-linearity of 
Eq. (p9|), as the difference ag — bg — 7. 

Fig. 12 The S^-weighted equivalent nucleon flux (j} eq {E^) is shown for different primary cosmic ray experiments 
[^8|"^ll . For each of these we plot the central value and the related error band. 

Fig. 13 The spectral index, 7 eo + 1, for the equivalent nucleon fluxes of Fig. 12, is shown for different primary cosmic 
ray experiments [p8|-^l|. For each of these we plot the central value and the related error band. 
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rric (GeV) 


£iF (TOT) 


{mr) 


11,1 )Vf D , , . 

a% NR (jib) 


°c* P (pb) 


1.1 


0.5 


2.53 


13.48 


13.5 ± 2.2 




1.0 


2.40 


13.48 






2.0 


2.10 


13.42 




1.2 


0.5 


1.46 


13.57 


13.5 ± 2.2 




1.0 


1.40 


13.54 






2.0 


1.23 


13.51 




1.25 


0.5 


1.18 


13.57 


13.5 ± 2.2 




1.0 


1.13 


13.54 






2.0 


1.00 


13.58 




1.3 


0.5 


0.96 


13.55 


13.5 ±2.2 




1.0 


0.92 


13.50 






2.0 


0.83 


13.53 




1.4 


0.5 


0.68 


13.51 


13.5 ±2.2 




1.0 


0.66 


13.51 






2.0 


0.61 


13.52 





TABLE I. Choice of parameters m c , /j,f and fin, that can reproduce the experimental total cross section a cB for charm 
production in pN collisions at 250 GeV from the E769 experiment. For each set of parameters, a™ N is the cross section 
calculated with the MNR program using MRST PDF. 
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Figure 2b. 
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Figure 12. 
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